function Step2_figure_4()
% *************************************************************************
% Simulations - for Figure 4
% cost shock in each country that generates a weighted average of 5% deflation in the country
% (its effect on other countries)

% 2008 linkage matrix used
% *************************************************************************
%%
clear; clc;
rootfolder      = input('Enter the location for the rootfolder in single quotes (''''): \n');
ppifolder       = [rootfolder '\analysis\PPI\MatlabFiles\'];
gammafolder     = [rootfolder '\analysis\WIOD\GAMMA\MatlabFiles\'];
toweightsfolder = [rootfolder '\analysis\WIOD\TOWEIGHTS\MatlabFiles\'];
totalOutputfolder = [rootfolder '\analysis\WIOD\totalOutput\MatlabFiles\'];
costshockfolder = [rootfolder '\analysis\COSTSHOCKS\17X31\'];
outputfolder = [rootfolder '\figures\simulations\'];
if exist(outputfolder,'dir') == 0
    mkdir(outputfolder);
end
cd(outputfolder);

%% Download Gamma matrix
GAMMAdata = importdata([gammafolder,'WIOD_gamma_2008.csv']);
Gamma_text = GAMMAdata.textdata(2:end,1:3);
GAMMA = GAMMAdata.data;

% List of countries & sectors
list_country = unique(Gamma_text(:,2),'stable');
n_country = size(list_country,1);
list_sector = unique(Gamma_text(:,3),'stable');
n_sector = size(list_sector,1);
n_cs = n_country * n_sector;

% Generate identity & 1 & 0 matrices
I = eye(size(GAMMA));
ONE = ones(size(GAMMA));
ZERO = zeros(size(GAMMA));
        
% Generate D matrix
D = diag(1 - sum(GAMMA,2));
D_inv = D \ I; 

% Generate  I-GAMMA and its inverse
IminusGamma = I - GAMMA ;
IminusGamma_inv = (I - GAMMA) \ I ;

% Download total output weights
TOWEIGHTSdata = importdata([toweightsfolder,'WIOD_TOweights_2008.csv']);
TOweights_text = TOWEIGHTSdata.textdata(2:end,1:3);
TOWEIGHTS = TOWEIGHTSdata.data;

% Download total output 
TOdata = importdata([totalOutputfolder,'WIOD_totalOutput_2008.csv']);
TO_text = TOdata.textdata(2:end,1:3);
totalOutput = TOdata.data;

% Check country - sector alignment across gamma & total output weights
gamma_sector = Gamma_text(:,3);
toweights_sector = TOweights_text(:,3);
assert(sum(strcmp(toweights_sector,gamma_sector)) == size(gamma_sector,1) );

nyears = 2011-1995+1;

% Cost shocks simulation
% III - countrefactuals: a cost shock in country iii that results in a 5%
% deflation in country iii -> what is the effect on other countries?

%create a selector array with selector for each country
selector = zeros(n_cs,n_country);
selector_not = zeros(n_cs,n_country);

%create a shock array
St_counterfactual = zeros(n_cs,n_country);
%create a corresponding dPPI array
dPPI_counterfactual = zeros(n_cs,n_country);
dPPI_counterfactual_direct = zeros(n_cs,n_country);

for iii = 1:n_country
    % select for each country, create shock, compute resulting dPPI
    temp_c_I = list_country(iii);
    selector(:,iii) = logical(strcmp(Gamma_text(:,2),temp_c_I) == 1);
    selector_not(:,iii) = logical(1-selector(:,iii));
    St_counterfactual(:,iii) = selector(:,iii) ; % .* (-0.05)
    dPPI_counterfactual(:,iii) = IminusGamma_inv * D * St_counterfactual(:,iii);

    % adjust because we want a 1% dPPI in the country %(-0.05)
    % compute the WEIGHTED mean of the shocked country
    domestic_weight = totalOutput(logical(selector(:,iii)));
    domestic_weight = domestic_weight ./ sum(domestic_weight);
    
    mean_temp = sum(domestic_weight .* dPPI_counterfactual(logical(selector(:,iii)),iii));
%     display('First run for country c_I =')
%     display(temp_c_I)
%     display('mean dppi for country')
%     display(mean_temp)
%     display('mean for other countries')
%     display( mean(dPPI_counterfactual(logical(selector_not(:,iii)),iii)) )
    
    % change shock accordingly
    St_counterfactual(:,iii) = selector(:,iii) .* abs(1/mean_temp) ;  % .* (-0.05) .* abs(0.05/mean_temp)
    dPPI_counterfactual(:,iii) = IminusGamma_inv * D *St_counterfactual(:,iii);
    
    domestic_weight = totalOutput(logical(selector(:,iii)));
    domestic_weight = domestic_weight ./ sum(domestic_weight);
    
    mean_temp = sum(domestic_weight .* dPPI_counterfactual(logical(selector(:,iii)),iii));
    
    display('Second run for country c_I =')
    display(temp_c_I)
    display('mean dppi for country')
    display(mean_temp)
    display('mean for other countries')
    display( mean(dPPI_counterfactual(logical(selector_not(:,iii)),iii)) )
    
    % now create the counterfactual dppi if there was only a linkage from
    % the source country to the destination country, an no indirect effect
    GAMMA_direct =  ZERO;
    GAMMA_direct(:,logical(selector(:,iii))) = GAMMA(:,logical(selector(:,iii)));
    
    D_direct = I - diag(sum(GAMMA_direct,2));
    IminusGamma_direct = I - GAMMA_direct;
    IminusGamma_inv_direct = (IminusGamma_direct*10000)\ I *10000;
    temp_direct = IminusGamma_inv_direct * D;
    
    dPPI_counterfactual_direct(:,iii) = temp_direct*St_counterfactual(:,iii);
end

% cost shock in the energy sector 10%
selector_energy = logical(strcmp(Gamma_text(:,3),'23') == 1);
selector_notEnergy = logical(1-selector_energy);
St_energy = selector_energy*10; % 0.1
dPPI_shock_energy = IminusGamma_inv * D *St_energy;
domestic_weight = totalOutput(logical(selector_energy));
domestic_weight = domestic_weight./ sum(domestic_weight);

mean_temp = sum(domestic_weight.*dPPI_shock_energy(logical(selector_energy)));
St_energy = selector_energy .* 10 .* abs(10/mean_temp); % .* (0.1) .* abs(0.1/mean_temp);
dPPI_shock_energy = IminusGamma_inv * D *St_energy;
domestic_weight = totalOutput(logical(selector_energy));
domestic_weight = domestic_weight./ sum(domestic_weight);
mean_temp = sum(domestic_weight.*dPPI_shock_energy(logical(selector_energy)));
% display('mean for energy')
% display(mean_temp)

%direct effect: shut down any linkage from non energy sector
GAMMA_direct =  zeros(size(GAMMA));
GAMMA_direct(:,logical(selector_energy)) = GAMMA(:,logical(selector_energy));
D_direct = eye(size(GAMMA_direct)) - diag(sum(GAMMA_direct,2));
IminusGamma_direct = eye(size(GAMMA_direct))-GAMMA_direct;
IminusGamma_inv_direct = (IminusGamma_direct*10000)\eye(size(GAMMA_direct))*10000;
temp_direct = IminusGamma_inv_direct*D;
dPPI_shock_energy_direct = temp_direct*St_energy;

% Export
numData = [totalOutput St_counterfactual dPPI_counterfactual dPPI_counterfactual_direct St_energy dPPI_shock_energy dPPI_shock_energy_direct];
numDataAsCell = num2cell(numData);

headers_St_counterfactual = strcat('cost_shock_',lower(list_country))';
headers_dPPI_counterfactual = strcat('Dppi_shock_comp_',lower(list_country))';
headers_dPPI_counterfactual_direct = strcat('Dppi_shock_dir_',lower(list_country))';
headers = ['total_output' headers_St_counterfactual headers_dPPI_counterfactual headers_dPPI_counterfactual_direct];
headers = [headers 'cost_shock_energy' 'Dppi_shock_comp_energy' 'Dppi_shock_dir_energy'];
headers = [GAMMAdata.textdata(1,1:size(Gamma_text,2)) headers];

outputData = [Gamma_text numDataAsCell];
output = [headers; outputData];

headerFmt = repmat('%s,',1,size(headers,2)-1);
textFmt = repmat('%s,',1,size(Gamma_text,2));
mathFmt = repmat('%1.6e,',1,size(numDataAsCell,2)-1);
fid = fopen(fullfile(outputfolder,strcat('outputMatlab_simulations.csv')), 'wt') ;
fprintf(fid,[headerFmt,'%s\n'], output{1,:});
for ll = 2:size(output,1)
    fprintf(fid,[textFmt,mathFmt,'%1.6e\n'], output{ll,:});
end
fclose(fid);

exit
end